Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :
La bioinfo c'est : <code>MERVEILLEUX</code>.
N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques/images pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.
Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :
mkdir -p ~/EvaluationM4M5/FASTQ
mkdir -p ~/EvaluationM4M5/CLEANING
mkdir -p ~/EvaluationM4M5/MAPPING
mkdir -p ~/EvaluationM4M5/QC
cd ~/EvaluationM4M5
tree ~/EvaluationM4M5
Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]
module load sra-tools
srun --cpus-per-task=6 fasterq-dump --split-files -p SRR10390685 --outdir FASTQ
Combien de reads sont présents dans les fichiers R1 et R2 ?
cd FASTQ
cat SRR10390685_1.fastq | echo $((`wc -l`/4))
cat SRR10390685_2.fastq | echo $((`wc -l`/4))
Les fichiers FASTQ contiennent 7066055(R1) + 7066055(R2) = 14132110 reads.
Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
cd ../MAPPING/
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz
Quelle est la taille de ce génome ?
zcat GCF_000009045.1_ASM904v1_genomic.fna.gz | awk '/^>/{if (l!="") print l; print; l=0; next}{l+=length($0)}END{print l}' |paste - -
La taille de ce génome est de 4215606 paires de bases.
Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz
Combien de gènes sont connus pour ce génome ?
zgrep -c $'\tgene\t' GCF_000009045.1_ASM904v1_genomic.gff.gz
4448 gènes sont recensés dans le fichier d’annotation.
Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit
cd ../
module load fastqc
gzip FASTQ/SRR10390685_1.fastq
gzip FASTQ/SRR10390685_2.fastq
srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_1.fastq.gz -o QC/ -t 8
srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_2.fastq.gz -o QC/ -t 8
La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?
car R2 comme le montre a plus de base dans la zone rouge que R1
Lien vers le rapport MulitQC
Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?
car il était nécessaire d’enlever “adapter content”
Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?
La profondeur de séquençage est de : 2.119142 G (total bases before filtering fastp) / 4215606 (ref genome size) = 502.7 (rather good) X.
Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp [3] qui vous semblent adéquats et justifiez-les.
module load fastp
srun --cpus-per-task 8 fastp --in1 FASTQ/SRR10390685_1.fastq.gz --in2 FASTQ/SRR10390685_2.fastq.gz --out1 CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz --out2 CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz --html CLEANING/fastp.html --thread 8 --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --json CLEANING/fastp.json
Les paramètres suivants ont été choisis :
| Parametre | Valeur | Explication |
|---|---|---|
Ces paramètres ont permis de conserver 13.554096M reads pairés, soit une perte de 95,90992427882319% des reads bruts.
Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [4] et samtools [5].
cd MAPPING/
GCF_000009045.1_ASM904v1_genomic.fna.gz > GCF_000009045.1_ASM904v1_genomic.fasta
module load bwa
srun bwa index GCF_000009045.1_ASM904v1_genomic.fasta
srun --cpus-per-task=32 bwa mem GCF_000009045.1_ASM904v1_genomic.fasta ../CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz ../CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz -t 32 > SRR10390685_on_GCF_000009045.1_ASM904v1.sam
module load samtools
srun --cpus-per-task=8 samtools view --threads 8 SRR10390685_on_GCF_000009045.1_ASM904v1.sam -b > SRR10390685_on_GCF_000009045.1_ASM904v1.bam
srun samtools sort SRR10390685_on_GCF_000009045.1_ASM904v1.bam -o SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam
srun samtools index SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam
srun samtools idxstats SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam > SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam.idxstats
srun samtools flagstat SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam > SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam.flagstat
samtools faidx GCF_000009045.1_ASM904v1_genomic.fasta
samtools faidx CP031214.1.fasta
Combien de reads ne sont pas mappés ?
13.57 M - 12.83M = 0.74 M reads ne sont pas mappés. Data from MulitQC report
13571369 - 12826829 = 744540 reads ne sont pas mappés. Data from flagstat and idxstats
Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [6]:
zgrep trmNF MAPPING/GCF_000009045.1_ASM904v1_genomic.gff.gz | awk '$3=="gene"' > trmNF.gff3
module load bedtools
bedtools bamtobed -i MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam > MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1.bed
bedtools intersect -a MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1.bed -b trmNF.bed -f 0.5 > result_trmNF_bed.txt
grep -c 'NC_000964.3' result_trmNF_bed.txt
2797 reads chevauchent le gène d’intérêt.
MulitQC report
export LC_ALL=en_US.UTF-8
export LANG=en_US.UTF-8
module load multiqc
srun multiqc -d . -o .
Utilisez IGV [7] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.
A work by Migale Bioinformatics Facility
https://migale.inrae.fr
Our two affiliations to cite us:
Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France
Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France